classdef QuesNewtonBasicMethod
    properties
         THRESHOLDGRADRMS = 10^-5;
         THRESHOLDMaxSTEP = 10^-5;
         THRESHOLDMaxGrad = 10^-5;
         THRESHOLDStepRMS = 10^-5;
         Rk = 0.2;        
         Grad
         eigVals
         eigVectors
         Dims
         MaxIterations = 1000;
         Scaler = 0.1;
         n_pos = 0;
         n_zero = 0; 
         n_neg = 0;
         E
         H 
         counter =0;
         NWMaxIterations = 10;         
         
    end
    methods
        function obj = QuesNewtonBasicMethod()
            obj.THRESHOLDGRADRMS = 10^-5;
            obj.THRESHOLDMaxSTEP = 10^-5;
            obj.THRESHOLDMaxGrad = 10^-5;
            obj.THRESHOLDStepRMS = 10^-5;
            obj.Rk = 0.2;
        end
        
        function flag = isATransitionState(obj, GradRMS, MaxGrad, ...
                MaxStep, StepRMS)
            flag  = false;
            if GradRMS <= obj.THRESHOLDGRADRMS &&  MaxGrad <= obj.THRESHOLDMaxSTEP && ...
                    MaxStep <= obj.THRESHOLDMaxGrad && StepRMS <=  obj.THRESHOLDStepRMS
                flag = true;
            end            
        end
        
        
        
        function Fit = EvaluateSchlegelLamba(obj,y)            
            Step = zeros(size(obj.Grad));
            for i = 1: obj.Dims                
                Grad_ = obj.eigVectors(:,i) * obj.Grad';
                Grad_ = Grad_/(obj.eigVals(i)- y ) ;                
                Step = Step - Grad_' * obj.eigVectors(:,i);
            end
            FinalFunc = 0;
            for i = 1: obj.Dims                
                FinalFunc = FinalFunc + Step(i)^2;
            end
            FinalFunc = sqrt(FinalFunc);
            Fit = (FinalFunc - obj.Rk)^2;
        end
        
    
        
        
        function [n_pos, n_zero, n_neg] = DetermineInertiaOfHessian(obj, ...
                eigVals, Dims)
            n_pos =0;
            n_zero = 0;
            n_neg = 0;
            for i=1: Dims
                if eigVals(i) == 0.0
                    n_zero = n_zero + 1;
                end
                if eigVals(i) > 0.0
                    n_pos = n_pos + 1;
                end
                if eigVals(i) < 0.0
                    n_neg = n_neg + 1;                    
                end                
            end
        end
        
        function Ce = CalculatePredictedChangeInEnergy(obj, Grad, Step, Hk)
            D = Step;
            G = Grad;
            Ce = G'*D + 0.5 * D' * Hk * D;
        end
        
        function [E, G, H] = CalculateRealEnergyDerivatives(obj, X)
            E = Function(X);
            G = Gradient(X);
            H = Hessian(X);
        end
        
        function Hk1 = UpdateHessian(obj, Step, Gk1, Gk, Hk)
            DT = Step';
            D = Step ;
            DTD = (DT * D);
            DDT = (D * DT);
            y = Gk1 - Gk - Hk * D;
            F = ((y' * D )*( DDT)) / DTD;
            Hk1 = Hk + (1/ DTD) * (y*DT + D*y' - F);
        end
        
        function [Step] = CalculateQausiStep(obj)
            Stp = QuasiNewtenStep;
            Step = Stp.CalculateStep(obj.Grad, obj.eigVals,obj.eigVectors,...
                obj.Dims);
        end
        
        function E = CalculateQuadraticEnergy(obj, Step, E1, Grad, H)
            E = E1 + Grad' * Step + 0.5 * Step' * H * Step;
        end
           
    end
end
